
require(emmeans)



##1 = Control, 2= Boy , 3= girl
model_TT <- lm(Therm_Trump ~ auth*as.factor(treatment), CFAdata_3)
anova(model_TT) #check the results
summary(model_TT$coefficients[2])
slopes_TT <- emtrends(model_TT, 'treatment', var = 'auth') #gets each slope
slopes_TT
pairs(slopes_TT)


model_BT <- lm(Therm_Biden ~ auth*as.factor(treatment), CFAdata_3)
anova(model_BT) #check the results

slopes_BT <- emtrends(model_BT, 'treatment', var = 'auth') #gets each slope
slopes_BT
pairs(slopes_BT)

model_DT <- lm(Therm_Dem ~ auth*as.factor(treatment), CFAdata_3)
anova(model_DT) #check the results

slopes_DT <- emtrends(model_DT, 'treatment', var = 'auth') #gets each slope
slopes_DT
pairs(slopes_DT)

model_GT <- lm(Therm_GOP ~ auth*as.factor(treatment), CFAdata_3)
anova(model_GT) #check the results

slopes_GT <- emtrends(model_GT, 'treatment', var = 'auth') #gets each slope
slopes_GT
pairs(slopes_GT)

model_BlmT <- lm(Therm_BLM ~ auth*as.factor(treatment), CFAdata_3)
anova(model_BlmT) #check the results

slopes_BlmT <- emtrends(model_BlmT, 'treatment', var = 'auth') #gets each slope

pairs(slopes_BlmT)

model_CT <- lm(Therm_Cops ~ auth*as.factor(treatment), CFAdata_3)
anova(model_CT) #check the results
summary(model_CT)
slopes_CT <- emtrends(model_CT, 'treatment', var = 'auth') #gets each slope
slopes_CT
pairs(slopes_CT)


##immigration
model_I <- lm(Immigration ~ auth*as.factor(treatment), CFAdata_3)
anova(model_I) #check the results

slopes_I <- emtrends(model_I, 'treatment', var = 'auth') #gets each slope
slopes_I
pairs(slopes_I)

model_G <- lm(GR ~ auth*as.factor(treatment), CFAdata_3)
anova(model_G) #check the results

slopes_G <- emtrends(model_G, 'treatment', var = 'auth') #gets each slope
slopes_G
pairs(slopes_G)

model_D <- lm(DP ~ auth*as.factor(treatment), CFAdata_3)
anova(model_D) #check the results

slopes_D <- emtrends(model_D, 'treatment', var = 'auth') #gets each slope
slopes_D
pairs(slopes_D)


model_HS <- lm(hostile_sexism ~ auth*as.factor(treatment), CFAdata_3)
anova(model_HS) #check the results
summary(model_HS)
slopes_HS <- emtrends(model_HS, 'treatment', var = 'auth') #gets each slope
slopes_HS
pairs(slopes_HS)



model_BS <- lm(ben_sexism ~ auth*as.factor(treatment), CFAdata_3)
anova(model_BS) #check the results

slopes_BS <- emtrends(model_BS, 'treatment', var = 'auth') #gets each slope
slopes_BS
pairs(slopes_BS)

model_HM <- lm(hostil_men ~ auth*as.factor(treatment), CFAdata_3)
anova(model_HM) #check the results

slopes_HM <- emtrends(model_HM, 'treatment', var = 'auth') #gets each slope
slopes_HM
pairs(slopes_HM)

model_BM <- lm(ben_men ~ auth*as.factor(treatment), CFAdata_3)
anova(model_BM) #check the results

slopes_BM <- emtrends(model_BM, 'treatment', var = 'auth') #gets each slope
slopes_BM
pairs(slopes_BM)

model_DA <- lm(auth ~ daughter*as.factor(treatment), CFAdata_3)
anova(model_DA) #check the results

slopes_DA <- emtrends(model_DA, 'treatment', var = 'daughter') #gets each slope
slopes_DA
pairs(slopes_DA)


model_PID <- lm(auth ~ pid*as.factor(treatment), CFAdata_3)
anova(model_PID) #check the results

slopes_PID <- emtrends(model_PID, 'treatment', var = 'pid') #gets each slope
slopes_PID
pairs(slopes_PID)

model_ID <- lm(auth ~ ideology*as.factor(treatment), CFAdata_3)
anova(model_ID) #check the results
summary(model_ID)
slopes_ID <- emtrends(model_ID, 'treatment', var = 'ideology') #gets each slope
slopes_ID
pairs(slopes_ID)






#convert the emtrends objects into something usable

TT <- summary(slopes_TT,infer=TRUE)
BT <- summary(slopes_BT,infer=TRUE)
DT <- summary(slopes_DT,infer=TRUE)
GT <- summary(slopes_GT,infer=TRUE)
BlmT <- summary(slopes_BlmT,infer=TRUE)
CT <- summary(slopes_CT,infer=TRUE)
Im <- summary(slopes_I,infer=TRUE)
Gr <- summary(slopes_G,infer=TRUE)
Dp <- summary(slopes_D,infer=TRUE)
Hs <- summary(slopes_HS,infer=TRUE)
Hm <- summary(slopes_HM,infer=TRUE)
Bs <- summary(slopes_BS,infer=TRUE)
Bm <- summary(slopes_BM,infer=TRUE)
Tp <- summary(slopes_PID,infer=TRUE)
Id <- summary(slopes_ID,infer=TRUE)
Da<- summary(slopes_DA,infer=TRUE)

Da$daughter.trend

model1FrameR <-data.frame(Variable = c(#"Daughter",
  "Partisanship",
  "Ideology"),
  
  Coefficient =c(
    Tp$pid.trend[1],
    Id$ideology.trend[1]#,
    #Da$daughter.trend[1]
  ),
  SE = c(
    Tp$SE[1],
    Id$SE[1]#,
    #Da$SE[1]
  ),
  modelName = c("Control"))

model2FrameR <-data.frame(Variable = c(#"Daughter",
  "Partisanship",
  "Ideology"),
  
  Coefficient =c(
    Tp$pid.trend[2],
    Id$ideology.trend[2]#,
    #Da$daughter.trend[2]
  ),
  SE = c(
    Tp$SE[2],
    Id$SE[2]#,
    #Da$SE[2]#
  ),
  modelName = c("Girl"))

model3FrameR <-data.frame(Variable = c(#"Daughter",
  "Partisanship",
  "Ideology"),
  
  Coefficient =c(
    Tp$pid.trend[3],
    Id$ideology.trend[3]#,
    #Da$daughter.trend[3]
  ),
  SE = c(
    Tp$SE[3],
    Id$SE[3]#,
    # Da$SE[3]
  ),
  modelName = c("Boy"))


# Combine these data.frames
allModelFrame <- data.frame(rbind(model1FrameR, 
                                  model2FrameR, model3FrameR))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Clean up the variable names

allModelFrame<-allModelFrame%>%
  mutate(Variable2 = recode_factor(.x=Variable,
                                   `(Intercept)` = "Intercept",
                                   `age`="Age",
                                   `bagain`="Born again",
                                   `ben_men`="Benevolent men",
                                   `ben_sexism`="Benevolent sexism",
                                   `chatt`="Church attendance",
                                   `education`="Education",
                                   `hostil_men`="Hostile men",
                                   `hostile_sexism`="Hostile sexism",
                                   `income2`="Income",
                                   `Male`="Male",
                                   `married`="Married",
                                   `noincome`="Income missing",
                                   `relimp`="Religion importance",
                                   `white`="White",
                                   `latino`="Latino",
                                   `daughter`="Daughter"
  ))%>%
  mutate(Model=modelName)
allModelFrame$Model<-allModelFrame$modelName
##drop the intercept.  Change the ordering of the variables
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
#zp1 <- zp1 + geom_linerange(aes(x = Variable2, ymin = Coefficient - SE*interval1,
#                                ymax = Coefficient + SE*interval2),
#                            lwd = 1/2, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/4, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()+labs(x="Variable", y="Coefficient")
zp1.PI <- zp1 + ggtitle("Correlates of authoritarianism by treatment\nParty and ideology")


print(zp1.PI)  # The trick to these is position_dodge().
predictors_fig<-zp1.PI
#thermometers

#convert the emtrends objects into something usable

TT <- summary(slopes_TT,infer=TRUE)
BT <- summary(slopes_BT,infer=TRUE)
DT <- summary(slopes_DT,infer=TRUE)
GT <- summary(slopes_GT,infer=TRUE)
BlmT <- summary(slopes_BlmT,infer=TRUE)
CT <- summary(slopes_CT,infer=TRUE)
Im <- summary(slopes_I,infer=TRUE)
Gr <- summary(slopes_G,infer=TRUE)
Dp <- summary(slopes_D,infer=TRUE)
Hs <- summary(slopes_HS,infer=TRUE)
Hm <- summary(slopes_HM,infer=TRUE)
Bs <- summary(slopes_BS,infer=TRUE)
Bm <- summary(slopes_BM,infer=TRUE)
Tp <- summary(slopes_PID,infer=TRUE)
Id <- summary(slopes_ID,infer=TRUE)




model1FrameR <-data.frame(Variable = c("Trump Therm", 
                                       "Biden Therm", 
                                       "Dem Therm",
                                       "GOP Therm",
                                       "BLM Therm",
                                       "Cops Therm"),
                          Coefficient =c(TT$auth.trend[1],
                                         BT$auth.trend[1],
                                         DT$auth.trend[1],
                                         GT$auth.trend[1],
                                         BlmT$auth.trend[1],
                                         CT$auth.trend[1]     ),
                          SE = c(TT$SE[1],
                                 BT$SE[1],
                                 DT$SE[1],
                                 GT$SE[1],
                                 BlmT$SE[1],
                                 CT$SE[1]),
                          modelName = c("Control"))

model2FrameR <-data.frame(Variable = c("Trump Therm", 
                                       "Biden Therm", 
                                       "Dem Therm",
                                       "GOP Therm",
                                       "BLM Therm",
                                       "Cops Therm"),
                          Coefficient =c(TT$auth.trend[2],
                                         BT$auth.trend[2],
                                         DT$auth.trend[2],
                                         GT$auth.trend[2],
                                         BlmT$auth.trend[2],
                                         CT$auth.trend[2]),
                          SE = c(TT$SE[2],
                                 BT$SE[2],
                                 DT$SE[2],
                                 GT$SE[2],
                                 BlmT$SE[2],
                                 CT$SE[2]),
                          modelName = c("Boy"))                         


model3FrameR <-data.frame(Variable = c("Trump Therm", 
                                       "Biden Therm", 
                                       "Dem Therm",
                                       "GOP Therm",
                                       "BLM Therm",
                                       "Cops Therm"),
                          Coefficient =c(TT$auth.trend[3],
                                         BT$auth.trend[3],
                                         DT$auth.trend[3],
                                         GT$auth.trend[3],
                                         BlmT$auth.trend[3],
                                         CT$auth.trend[3]),
                          SE = c(TT$SE[3],
                                 BT$SE[3],
                                 DT$SE[3],
                                 GT$SE[3],
                                 BlmT$SE[3],
                                 CT$SE[3]),
                          modelName = c("Girl"))
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1FrameR, 
                                  model2FrameR, model3FrameR))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Clean up the variable names

allModelFrame<-allModelFrame%>%
  mutate(Variable2 = recode_factor(.x=Variable,
                                   `(Intercept)` = "Intercept",
                                   `age`="Age",
                                   `bagain`="Born again",
                                   `ben_men`="Benevolent men",
                                   `ben_sexism`="Benevolent sexism",
                                   `chatt`="Church attendance",
                                   `education`="Education",
                                   `hostil_men`="Hostile men",
                                   `hostile_sexism`="Hostile sexism",
                                   `income2`="Income",
                                   `Male`="Male",
                                   `married`="Married",
                                   `noincome`="Income missing",
                                   `relimp`="Religion importance",
                                   `white`="White",
                                   `latino`="Latino",
                                   `daughter`="Daughter"
  ))%>%
  mutate(Model=modelName)
allModelFrame$Model<-allModelFrame$modelName
##drop the intercept.  Change the ordering of the variables
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
#zp1 <- zp1 + geom_linerange(aes(x = Variable2, ymin = Coefficient - SE*interval1,
#                                ymax = Coefficient + SE*interval2),
#                            lwd = 1/2, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/4, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()+labs(x="Variable", y="Coefficient")
zp1 <- zp1 + ggtitle("Correlates of authoritarianism by treatment\nFeeling thermometers")


therm_fig<-print(zp1)  # The trick to these is position_dodge().
therm_fig



##ISSUES



model1FrameR <-data.frame(Variable = c("Immigration",
                                       "Gay Rights",
                                       "Death Penalty"),
                          Coefficient =c(Im$auth.trend[1],
                                         Gr$auth.trend[1],
                                         Dp$auth.trend[1]),
                          SE = c(Im$SE[1],
                                 Gr$SE[1],
                                 Dp$SE[1]
                          ),
                          modelName = c("Control"))

model2FrameR <-data.frame(Variable = c("Immigration",
                                       "Gay Rights",
                                       "Death Penalty"),
                          Coefficient =c(Im$auth.trend[2],
                                         Gr$auth.trend[2],
                                         Dp$auth.trend[2]
                          ),
                          SE = c(
                            Im$SE[2],
                            Gr$SE[2],
                            Dp$SE[2]),
                          modelName = c("Boy"))                         


model3FrameR <-data.frame(Variable = c("Immigration",
                                       "Gay Rights",
                                       "Death Penalty"),
                          Coefficient =c(
                            Im$auth.trend[3],
                            Gr$auth.trend[3],
                            Dp$auth.trend[3]),
                          SE = c(
                            Im$SE[3],
                            Gr$SE[3],
                            Dp$SE[3]),
                          modelName = c("Girl"))
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1FrameR, 
                                  model2FrameR, model3FrameR))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Clean up the variable names

allModelFrame<-allModelFrame%>%
  mutate(Variable2 = recode_factor(.x=Variable,
                                   `(Intercept)` = "Intercept",
                                   `age`="Age",
                                   `bagain`="Born again",
                                   `ben_men`="Benevolent men",
                                   `ben_sexism`="Benevolent sexism",
                                   `chatt`="Church attendance",
                                   `education`="Education",
                                   `hostil_men`="Hostile men",
                                   `hostile_sexism`="Hostile sexism",
                                   `income2`="Income",
                                   `Male`="Male",
                                   `married`="Married",
                                   `noincome`="Income missing",
                                   `relimp`="Religion importance",
                                   `white`="White",
                                   `latino`="Latino",
                                   `daughter`="Daughter"
  ))%>%
  mutate(Model=modelName)
allModelFrame$Model<-allModelFrame$modelName
##drop the intercept.  Change the ordering of the variables
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
#zp1 <- zp1 + geom_linerange(aes(x = Variable2, ymin = Coefficient - SE*interval1,
#                                ymax = Coefficient + SE*interval2),
#                            lwd = 1/2, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/4, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()+labs(x="Variable", y="Coefficient")
zp1 <- zp1 + ggtitle("Correlates of authoritarianism by treatment\nIssue positions")


print(zp1)  # The trick to these is position_dodge().
issue_fig<-zp1


##SEXISM


model1FrameR <-data.frame(Variable = c("Hostile Sexism",
                                       "Hostile Men",
                                       "Benevolent Sexism",
                                       "Benevolent Men"
),
Coefficient =c(
  Hs$auth.trend[1],
  Hm$auth.trend[1],
  Bs$auth.trend[1],
  Bm$auth.trend[1]),
SE = c(
  Hs$SE[1],
  Hm$SE[1],
  Bs$SE[1],
  Bm$SE[1]),
modelName = c("Control"))

model2FrameR <-data.frame(Variable = c(
  "Hostile Sexism",
  "Hostile Men",
  "Benevolent Sexism",
  "Benevolent Men"),
  Coefficient =c(Hs$auth.trend[2],
                 Hm$auth.trend[2],
                 Bs$auth.trend[2],
                 Bm$auth.trend[2]),
  SE = c(Hs$SE[2],
         Hm$SE[2],
         Bs$SE[2],
         Bm$SE[2]),
  modelName = c("Boy"))                         


model3FrameR <-data.frame(Variable = c("Hostile Sexism",
                                       "Hostile Men",
                                       "Benevolent Sexism",
                                       "Benevolent Men"),
                          Coefficient =c(Hs$auth.trend[3],
                                         Hm$auth.trend[3],
                                         Bs$auth.trend[3],
                                         Bm$auth.trend[3]),
                          SE = c(Hs$SE[3],
                                 Hm$SE[3],
                                 Bs$SE[3],
                                 Bm$SE[3]),
                          modelName = c("Girl"))
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1FrameR, 
                                  model2FrameR, model3FrameR))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Clean up the variable names

allModelFrame<-allModelFrame%>%
  mutate(Variable2 = recode_factor(.x=Variable,
                                   `(Intercept)` = "Intercept",
                                   `age`="Age",
                                   `bagain`="Born again",
                                   `ben_men`="Benevolent men",
                                   `ben_sexism`="Benevolent sexism",
                                   `chatt`="Church attendance",
                                   `education`="Education",
                                   `hostil_men`="Hostile men",
                                   `hostile_sexism`="Hostile sexism",
                                   `income2`="Income",
                                   `Male`="Male",
                                   `married`="Married",
                                   `noincome`="Income missing",
                                   `relimp`="Religion importance",
                                   `white`="White",
                                   `latino`="Latino",
                                   `daughter`="Daughter"
  ))%>%
  mutate(Model=modelName)
allModelFrame$Model<-allModelFrame$modelName
##drop the intercept.  Change the ordering of the variables
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
#zp1 <- zp1 + geom_linerange(aes(x = Variable2, ymin = Coefficient - SE*interval1,
#                                ymax = Coefficient + SE*interval2),
#                            lwd = 1/2, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/4, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()+labs(x="Variable", y="Coefficient")
zp1 <- zp1 + ggtitle("Correlates of authoritarianism by treatment\nSexism measures")


print(zp1)  # The trick to these is position_dodge().
sexism_fig<-zp1


Tp <- summary(slopes_PID,infer=TRUE)
Id <- summary(slopes_ID,infer=TRUE)
Da<- summary(slopes_DA,infer=TRUE)

Da$daughter.trend

model1FrameR <-data.frame(Variable = c("Daughter",
                                       "Partisanship",
                                       "Ideology"),
                          
                          Coefficient =c(
                            Tp$pid.trend[1],
                            Id$ideology.trend[1],
                            Da$daughter.trend[1]),
                          SE = c(
                            Tp$SE[1],
                            Id$SE[1],
                            Da$SE[1]),
                          modelName = c("Control"))

model2FrameR <-data.frame(Variable = c("Daughter",
                                       "Partisanship",
                                       "Ideology"),
                          
                          Coefficient =c(
                            Tp$pid.trend[2],
                            Id$ideology.trend[2],
                            Da$daughter.trend[2]),
                          SE = c(
                            Tp$SE[2],
                            Id$SE[2],
                            Da$SE[2]),
                          modelName = c("Boy"))

model3FrameR <-data.frame(Variable = c("Daughter",
                                       "Partisanship",
                                       "Ideology"),
                          
                          Coefficient =c(
                            Tp$pid.trend[3],
                            Id$ideology.trend[3],
                            Da$daughter.trend[3]),
                          SE = c(
                            Tp$SE[3],
                            Id$SE[3],
                            Da$SE[3]),
                          modelName = c("Girl"))


# Combine these data.frames
allModelFrame <- data.frame(rbind(model1FrameR, 
                                  model2FrameR, model3FrameR))  # etc.

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

#Clean up the variable names

allModelFrame<-allModelFrame%>%
  mutate(Variable2 = recode_factor(.x=Variable,
                                   `(Intercept)` = "Intercept",
                                   `age`="Age",
                                   `bagain`="Born again",
                                   `ben_men`="Benevolent men",
                                   `ben_sexism`="Benevolent sexism",
                                   `chatt`="Church attendance",
                                   `education`="Education",
                                   `hostil_men`="Hostile men",
                                   `hostile_sexism`="Hostile sexism",
                                   `income2`="Income",
                                   `Male`="Male",
                                   `married`="Married",
                                   `noincome`="Income missing",
                                   `relimp`="Religion importance",
                                   `white`="White",
                                   `latino`="Latino",
                                   `daughter`="Daughter"
  ))%>%
  mutate(Model=modelName)
allModelFrame$Model<-allModelFrame$modelName
##drop the intercept.  Change the ordering of the variables
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
#zp1 <- zp1 + geom_linerange(aes(x = Variable2, ymin = Coefficient - SE*interval1,
#                                ymax = Coefficient + SE*interval2),
#                            lwd = 1/2, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/4, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()+labs(x="Variable", y="Coefficient")
zp1 <- zp1 + ggtitle("Correlates of authoritarianism by treatment\nPredictors")


print(zp1)  # The trick to these is position_dodge().

library(cowplot)
pdf(file="Figure9.pdf")

plot_grid(predictors_fig ,therm_fig, issue_fig  )
dev.off()
